Cite this Article
Liu Nan, Wang Dan-Ni, Liu Hai-Ying, Yang Hong-Li, Yang Lian-Gui. Potential dynamic analysis of tumor suppressor p53 regulated by Wip1 protein. Chinese Physics B, 2020, 29(6): 068704  
Permissions
Potential dynamic analysis of tumor suppressor p53 regulated by Wip1 protein
Liu Nan, Wang Dan-Ni, Liu Hai-Ying, Yang Hong-Li †, Yang Lian-Gui
School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China

 

† Corresponding author. E-mail: imuyhl@imu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 11762011).

Abstract

The tumor suppressor p53 plays a key role in protecting genetic integrity. Its dynamics have important physiological significance, which may be related to the cell fate. Previous experiments have shown that the wild-type p53-induced phosphatase 1 (Wip1) protein could maintain p53 oscillation. Therefore, we add Wip1 to remodel the p53 network. Firstly, we use the binomial τ-leap algorithm to prove our model stable under internal noise. Then, we make a series of bifurcation diagrams, that is, p53 levels as a function of p53 degradation rate at different Wip1 generation rates. The results illustrate that Wip1 is essential for p53 oscillation. Finally, a two-dimensional bifurcation diagram is made and the stability of some p53 dynamics under external noise is analyzed by potential landscape. Our results may have some implications for artificially interfering with p53 dynamics to achieve tumor suppression.

PACS: ;87.85.Xd;;87.57.cm;;82.40.Bj;;82.39.Rt;
Keyword:;Wip1;;p53-dynamics;;oscillation;
1. Introduction

The transcription factor p53 plays an important role in cell growth, senescence, DNA repair, and programmed cell death.[1–4] Exploring the dynamics of p53 may be helpful for understanding its fine-tuning of cellular process. Researchers have proposed that p53 performs its functions in relation to its dynamic behavior.[5] For example, transient p53 pulses promote DNA repair, while jumping to a highly steady state may cause apoptosis,[6] or apoptosis requires continuous p53 pulses.[7,8]

Wild-type p53-induced phosphatase 1 (Wip1) is one of the oncoproteins that inhibit p53 function. It have been found that Wip1 has overexpression in multiple human cancers such as breast, lung, pancreas, bladder, liver cancer, and in neuroblastomas.[9–12] In contrast, Wip1 gene knockout mice are resistant to oncogene-induced cancers,[11] suggesting that inhibition of Wip1 activity may be beneficial for cancer treatment. It is still necessary to further explore how Wip1 mediates the p53 function and the p53 dynamics. In this paper, we investigate how Wip1 affects the p53 pulse dynamic under cellular DNA double strand breaks (DSB).

The oscillating behavior of p53 has attracted attention after the experiment of Lev et al.[13] They used immunoblotting to detect p53 damped oscillations in population cells under ionizing radiation. Later, Lahav et al.[14] observed fluorescently labeled p53 protein in single cells, and found that p53 in single cells exhibits long-lasting undamped oscillations. Batchelor et al.[15] confirmed that Wip1 is the key to the sustained pulse behavior of p53. Zhang et al.[16] and Sun et al.[17] added p53-Wip1-ATM (the ataxia mutant) negative feedback based on the p53-Mdm2 model, indicating that the p53 pulse is driven by the upstream ATM pulse.

In this article, we rebuilt a p53 model that contains only p53, Mdm2, ATM and Wip1, and prove the correctness of our model by the time occurring diagram. Then we use the binomial τ-leap algorithm in the simulation. The results show that our model is quite robust under internal noise. Using the bifurcation diagram can effectively analyze the potential dynamics of the system.[18] We fix the maximum generation rate of Wip1 and plot a series of bifurcation diagrams of p53 steady level vs the maximum degradation rate of p53, and find a variety of dynamics in the p53 network, including monostable, oscillation, coexistence of monostable and oscillation, and bistable behavior, etc.

Finally, we make a two-dimensional bifurcation diagram, and find the parameter area corresponding to each kinetic. The potential landscape theory is used to study the stochastic properties of biological signal networks under external noise,[19–21] which demonstrate that there exists the p53 potential dynamics mode when the parameters are determined. The landscape suggests that each dynamics mode of p53 is quite stable under fit external noise.

2. Model

It is easy to find that p53-Mdm2 negative feedback loop (NFL) can generate oscillations by supplying time delay to Mdm2 synthesis term.[22,23] Adding positive feedback loop (PFL) to the p53-Mdm2 NFL can lead to supercritical Hopf bifurcation, which is more in line with the response manner of p53 ‘all or none’ under pressure.[24] Although many models have explained the effect of Wip1 on p53 oscillation, few of them consider the effect of the combined Wip1 and p53 degradation rate on the p53 kinetics. For this reason, we built a simplified model, and use the hidden time delay method, i.e., divide Mdm2 into nucleus and cytoplasm. We also introduce phosphorylated Mdm2 to promote the formation of p53, so there are one PFL and two NFLs in our model, which is similar to the model of Zhuge et al.[25]

In most cases, p53 is suppressed to a low level by Mdm2 in the nucleus.[26] Our model is based on the case of DNA continuous damage. Under this condition, ATM sensor can always sense the damage signal and be activated (phosphorylated) continuously.[27] Activated ATM further promotes the phosphorylation of p53 and Mdm2.[28,29] The phosphorylation of p53 is conducive to its stable accumulation, while the phosphorylation of Mdm2 accelerates the self-ubiquitination degradation and promotes the synthesis of p53.[29,30]

Damage signals can help p53 accumulation, but p53 will inhibit the increase of itself through NFL. On the one hand, p53 promotes Mdm2 transcription to form p53-Mdm2 NFL. On the other hand, it promotes Wip1 transcription to form p53-Wip1-ATM NFL. However, the existence of phosphorylated form of Mdm2 makes not only NFL but also PFL in p53-Mdm2. The model frame is shown in Fig. 1, which contains four proteins, p53, ATM, Wip1, and Mdm2. We divide Mdm2 into three forms, i.e., nuclear form (Mdm2n), cytoplasmic form (Mdm2c), and phosphorylated form (Mdm2p).

Fig. 1. Schematic diagram of the network model. Promotion and state transition is represented by arrow-headed lines, inhibition is denoted by bar-headed lines. Here 0,…,6 represent seven reaction channels using the Hill function.

All promotion and inhibition are represented by the Hill function. Here ρi (i = 0, …, 6) represent the maximum generation (phosphorylation) or degradation (dephosphorylation) proportional constant. We set ρi (i = 0, …, 6) ≥ 90 % to ensure that each component is mainly controlled by members within the network; ki (i = 0, …, 6) are the half-saturated concentration, where the subscript i represent the reaction channel, as shown in Fig. 1. Since the true concentration unit cannot be estimated by the experimental data, the concentration of all proteins in the model are dimensionless. By comparing the oscillation period with the experiment, the unit of time can be given in minutes. The corresponding ordinary differential equations (ODEs) of the model are given in the following, in which [⋅] denotes concentration of components.

ATMp reads

It contains activation and inactivation terms. Here the inactivation term is Wip1-dependent.

For the equation of p53, its generation and activation are controlled by Mdm2p and ATMp, respectively, so we multiply the two Hill functions as the active p53 increase term, and the p53 decrease term is only related to Mdm2n,

For the equation of Mdm2 (Eqs. (3)–(5)), p53 promotes the production of Mdm2c, and ATMp promotes the conversion of Mdm2c to Mdm2p. These two effects are expressed by Hill functions. Mdm2 phosphorylation (dephosphorylation) rate is kp (kq), and its rate of entry (exit) to the nucleus is marked as kin (kout). The degradation rate of Mdm2 depends on the form of Mdm2. The degradation rate of Mdm2p is greater than Mdm2c and assuming g0 times that of Mdm2c. For simplicity, our model does not include the phosphorylated form of Mdm2 in the nucleus, and assumes that the degradation rate of Mdm2n is linearly related to ATMp with a proportional constant r,

Wip1’s equation (6) includes a p53-dependent generation term and a steady-rate degradation term,

Here v and d represent the maximum synthesis (activation) and degradation (inactivation) rates, and the subscripts are the protein species. Most of the parameters are from references, a small number of parameters are slightly modified in order to reproduce the dynamic behavior of p53. We choose all time-varying protein concentrations to be 0 as the initial conditions. Bifurcation analysis was implemented using the free software Xppaut, and other numerical analysis was performed using software Matlab. The parameters are listed in Table 1.

Table 1.

Simulation parameters.

.
3. Results
3.1. System oscillation in the ODEs’ solution and τ-leap simulation

The numerical solution of the ODEs is shown by the solid line in Fig. 2. The period of the p53 oscillation is about 5.5 hours, which is consistent with the experimental data.[14] Oscillation behavior is not unique to p53 and Mdm2, but is present throughout the entire network of proteins, including ATM and Wip1. Our model suggest that the oscillation of p53 is driven by the oscillation of ATM (the solid line in Fig. 2), and Wip1 works to maintain the pulse dynamics (compare the solid line and the dashed line in Fig. 2), consistent with the experiment.[15] The importance of Wip1 to system oscillation in this model will be analyzed and discussed delicately in the later subsection.

Fig. 2. Time courses of protein concentrations with Wip1 (solid line), without Wip1 (dashed line). [Mdm2] represents the total concentration of three forms of Mdm2.

If the number of molecules of a protein is not very large, or the concentration of this protein is very low, ODEs are no longer applicable to the corresponding model, and random algorithms that characterize internal noise should be used. Stochastic simulations were performed using the binomial τ-leap algorithm instead of the Gillespie algorithm. That is, biochemical reactions are divided into two types, for the first type of reaction (only involving protein production), the number of times that each reaction occurs follows a Poisson distribution, and the average value is equal to the firing rate multiplied by the time step, for another kind of reaction (involving protein conversion or degradation), replace the Poisson distribution with a Binomial distribution (see Ref. [31] for details).

The binomial τ-leap algorithm is similar to the deterministic model when the cell volume Ω is large enough (compare Figs. 3(a)–3(c)). Our model under binomial τ-leap simulation will have recognizable pulse dynamics even if the cell volume is very small (for example Ω = 103). Therefore, it can be seen that the model is robust and reliable.

Fig. 3. Time courses of [p53]. Dimensionless cell volume Ω is (a) 103, (b) 104, and (c) 105.
3.2. Possible dynamics of the system at several levels Wip1 regulation

In the following we focus on stable steady states and limit cycles because these solutions shape time-dependent responses. As we all know, the p53-Mdm2 NFL is the basis of p53 oscillation, and the strength of this NFL is directly controlled by the degradation parameter dp53. Bifurcation analysis in Fig. 4 reveals the dependence of steady-state p53 levels on the dp53. In this section, we consider four values of vwip1.

Fig. 4. Bifurcation diagram of p53 vs p53 degradation rate (dp53): (a) Wip1 synthesis rate vwip1 = 0, (b) vwip1 = 0.046, (c) vwip1 = 0.051 and (d) vwip1 = 0.12. The stable and unstable steady states are indicated by red and black lines, respectively. Green circles and blue circles are the maxima and minima of the stable and unstable limit cycles, respectively.

When vwip1 = 0, there are two saddle node (SN) points in Fig. 4(a). Conversion between the low and high expressions of p53 is realized via bistable switching. Only when dp53 exceeds the critical value at SN2 can p53 switch to low levels. Conversely, when dp53 reduces to the value at SN1, p53 can switch to high levels. The system exhibits typical bistability (two stable steady states and one unstable steady state) for dp53 ∈ (SN1, SN2). This switch phenomenon is similar to Ref.[16].

We fix the value of vwip1 = 0.046 (Fig. 4(b)), for dp53 ∈ (0, SN1) or dp53 ∈ (SN2, + ∞), the system exhibits monostability (one stable steady state). Between SN1 and SN2, the subcritical hopf bifurcation point (HB), homoclinic bifurcation point (HC), and limit cycle fold bifurcations point (LPC) appear, the system exhibits typical bistability for dp53 ∈ (SN1, LPC). Atypical tristability (the stable limit cycle coexists with two stable steady state, one unstable limit cycle and one unstable steady states) occurs for dp53 ∈ (LPC, HC). For dp53 ∈ (HC, HB), the system returns to typical bistability (two stable steady states coexists with one unstable steady state and one unstable limit cycle). When the value of dp53 arises at HB, upper branch becomes unstable and system emerges excitability (one stable steady states and two unstable steady states) for dp53 ∈ (HB, SN2).

The character of the bifurcation diagram for vwip1 = 0.051 (Fig. 4(c)) is qualitatively similar to that for vwip1 = 0.046 (Fig. 4(b)). However, it is visually clearer that the distinguish between Fig. 4(c) and Fig. 4(b) is whether the HC exists. When vwip1 = 0.051, HC disappears and SN1 changes into a saddle note invariant circle bifurcation point (SNIC), in this case, atypical bistability (the stable limit cycle coexists with one stable steady state and one unstable limit cycle) appears for dp53 ∈ (LPC, SNIC), and typical bistability appears for dp53 ∈ (SNIC, HB).

As shown in Fig. 4(d), with dp53 increasing (for fixed vwip1 = 0.12) the system proceeds sequentially through five different dynamics in the bifurcation diagram. First of all, experience high steady state. The oscillation region is constrained by LPC and SNIC, between HB and SNIC, the system into the oscillatory state (stable limit cycle and unstable steady state), and then enters the excited region, finally the system stays at the low steady state. The oscillation appears in large amplitude, similar to Ref.[24]. Clearly, the long-term dynamic behaviors of [p53] in the above figures well coincide with these features.

3.3. Codimension-two bifurcation analysis

In Fig. 5 we display the two-dimensional bifurcation diagram in the (dp53, vwip1)-plane. As shown, the system involves five types of bifurcations: HB, SN, SNIC, LPC, and HC. HC bifurcations line arises at one SNIC point (dp53, vwip1) = (0.7727, 0.05).

Fig. 5. Two-dimensional bifurcation diagram showing various types of bifurcation lines (dp53, vwip1)-plane. The bifurcation lines divide the parameter domain into seven subdomains R1, …, R7.

The bifurcation lines divide the (dp53, vwip1)-plane into seven subregions in which the system exhibits different behaviors:

monostability (one stable steady state).

atypical bistability (one stable steady state coexists with stable limit cycle and unstable limit cycle).

oscillatory (stable limit cycle and one unstable steady state).

monostability/excitebility (one stable steady state coexists with two unstable steady states).

typical bistability (two stable steady state coexists with unstable limit cycle and unstable steady state).

typical bistability (two stable steady states and one unstable steady state).

atypical tristability (two stable steady states coexists with stable limit cycle and unstable limit cycle).

Our model omits the upstream and downstream modules of p53. However, with the help of other models, it can be known that the oscillation or multistable state of p53 is the prerequisite for a correct response of cells under pressure. For example, in the “digital” cell fate determination model, the number of p53 pulses corresponds to the degree of cell damage, and irreparable cells are removed in time in this mode.[16] However, in the “analog” cell fate decision model, the multi steady stable state of p53 is indispensable for cell fate determination.[32] A moderate level of p53 is conducive to DNA repair and a high level is conducive to apoptosis. In summary, cells need to avoid the R1 region to make a proper pressure response.

3.4. Potential landscape analysis

The dynamics of biological signal systems can be described stochastically, dx/dt = F(x) + ζ. Here F indicates the driving force describing the dynamics of the system, the form of F(x) is displayed in the Eqs. (1)–(6). Here x indicates the concentration of different components in the system, ζ is external noise term, its statistical properties can usually be assumed Gaussian and white, i.e., 〈ζ(t)ζ(t′)〈 = 2Dδ (t – t′), where δ (t) is the Dirac function. D is the matrix of the diffusion coefficient. It is assumed that D is homogeneous and constant. The probabilistic development of the corresponding dynamic system can be described by the diffusion equation in the form , where represents probability flux.[19]

It is hard to solve the high dimensional diffusion equations directly, we give a uniform distribution of the initial probability P(x,0) using statistical methods to evolve the P(x,t) over time. When the distribution of P becomes stable with time t increasing, one can get steady state of probability distribution, Pss. Define potential function U = –ln P,[19] then we can obtain the image of the potential landscape at steady state, Uss, as shown in Fig. 6.

Fig. 6. Three-dimensional potential landscape. The parameters are taken from the parameter regions: (a) R1, (b) R3, (c) R2, (d) R6. The diffusion coefficient D is 10–6. The height of potential (Uss) see the color bar on the right.

When the parameters are in the R1 region, the potential landscape is funnel-shaped (Fig. 6(a)), and the global minimum corresponds to the steady-state solution of the system. When the parameter is taken from the R3 region, a closed loop valley appears on the potential landscape (Fig. 6(b)), and the system is attracted to the closed loop instead of one point. The depth of the valley is inhomogeneous, because the system passes through each state at a different rate. When the system coexists with steady state and oscillation, both the funnel and the valley appear on the landscape (Fig. 6(c)). When the system is bistable, two funnel shapes appear (Fig. 6(d)), and the two local minima correspond to the two steady states of the system. Therefore, each p53 kinetic mode controlled by Wip1 expression rates is rather stable at D = 10–6.

Once the noise intensity D becomes large, the closed loop valley in oscillation mode turns into a basin (compare Figs. 7(a)–7(c)), indicating that the closed-loop attraction is weak, the system is easy to escape the oscillation state. In short, p53 oscillations may be stable under external noise in a real cellular environment.

Fig. 7. The two-dimensional landscape. The diffusion coefficient D is (a) 10–6, (b) 10–4, (c) 10–3. For the height of potential (Uss) see the color bar on the right.
4. Conclusion

We have established a model of the p53 network system and analyzed the dynamic mechanism. Using bifurcation analysis, we find that sufficient Wip1 generation rate and appropriate p53 degradation rate can cause the p53 system to oscillate. Wip1 is of great significance for maintaining p53 oscillation, this is consistent with the phenomenon observed in the experiment.[15] We further use codimension-two bifurcation analysis to explore the combined effect of the two parameters on the potential kinetic mechanism of the p53 system, and find different dynamic behaviors in different regions of the parametric plane. The important influence of Wip1 on the dynamic diversity of p53 network is verified.

Cells respond appropriately to pressure, both vwip1 and dp53 need to take appropriate values so that p53 avoids monostable regions. In the p53 monostable region, low steady state can cause cells to fail to apoptosis, even in severe damage. High steady state may make the cell sensitive excessively to the damage signal, causing the cell to die during physiological DNA damage fluctuations (such as mitosis). We can note that vwip1 is too large, leading to the parameters to be in the R1 region on the upper right of Fig. 5, which may be one reason for cell canceration. Therefore, the possibility of clinical treatment of cancer can be considered from the aspect of reducing vwip1 or dp53.

We use the potential landscape and observe the distribution of different protein states under several parameter values, which is consistent with the bifurcation diagram. The robustness of p53 oscillations under external noise is explained by the landscape. Since the dynamics of p53 is inseparable from cell fate, we unveiled the p53 rich dynamic mechanism, hoping to contribute to cancer treatment, and we expect to utilize mathematical models in this way to build more accurate and predictive simulators.

Reference
[1] Vousden K H Prives C 2009 Cell 137 413
[2] Kastan M B Onyekwere O Sidransky D Vogelstein B Craig R W 1991 Cancer Res. 51 6304
[3] Sengupta S Harris C C 2005 Nat. Rev. Mol. Cell Biol. 6 44
[4] Elmore S 2007 J. Toxicol Pathol 35 495
[5] Purvis J E Karhohs K W Mock C Batchelor E Loewer A Lahav G 2012 Science 336 1440
[6] Zhang Q H Tian X J Liu F Wang W 2014 FEBS Lett. 588 4369
[7] Zhang X P Liu F Cheng Z Wang W 2009 Proc. Natl. Acad. Sci. USA 106 12245
[8] Sun T Chen C Wu Y Zhang S Cui J Shen P 2009 BMC Bioinformatics 10 190
[9] Li J Yang Y Peng Y et al. 2002 Nat. Genetics 31 133
[10] Rauta J Alarmo E L Kauraniemi P Karhu R Kuukasjärvi T Kallioniemi A 2006 Breast Cancer Res. Treat. 95 257
[11] Hirrison M Li J Degenhardt Y Hoey T Powers S 2004 Trends Mol. Med. 10 359
[12] Bulavin D V Demidov O N Saito S et al. 2002 Nat. Genet. 31 210
[13] Lev B R Maya R Segel L A Alon U Levine A J Oren M 2000 Proc. Natl. Acad. Sci. USA 97 11250
[14] Lahav G Rosenfeld N Sigal A Geva-Zatorsky N Levine A J Elowitz M B Alon U 2004 Nat. Genet. 36 147
[15] Batchelor E Mock C S Bhan I Loewer A Lahav G 2008 Mol. Cell 30 277
[16] Zhang X P Liu F Wang W 2011 Proc. Natl. Acad. Sci. USA 108 8990
[17] Sun T Z Yang W Liu J Shen P P 2011 PLoS. ONE 6 e27882
[18] Wang D G Zhou C H Zhang X P 2017 Chin. Phys. 26 128709
[19] Wang J Li X Wang E 2008 Proc. Natl. Acad. Sci. USA 105 12271
[20] Li C Ye L 2019 J. Chem. Phys. 151 175101
[21] Bi Y H Yang Z Q He X Y 2016 Acta. Phys. Sin. 65 028701 in Chinese
[22] Xia J F Jia Y 2010 Chin. Phys. 19 040506
[23] Zhuo Y Z Yan S W Zhang L J 2007 Acta. Phys. Sin. 56 2442 in Chinese
[24] Zhang T Brazhnik P Tyson J J 2007 Cell Cycle 6 85
[25] Zhuge C Sun X Chen Y Lei J 2016 J. Theor. Biol. 388 1
[26] Levine A J 1997 Cell 88 323
[27] Falck J Coates J Jackson S P 2005 Nature 434 605
[28] Prives C 1998 Cell 95 5
[29] Stommel J M Wahl G M 2004 EMBO J. 23 1547
[30] Hamard P J Manfredi J J 2012 Cancer Cell 21 3
[31] Wang D G Wang S B Huang B Liu F 2019 Sci. Rep. 9 5883
[32] Zhang X P Liu F Wang W 2012 Biophys. J. 102 2251